```{r}

library(pacman)
p_load("fixest", "haven", "conleyreg", "readxl", "dplyr", "lmtest", "ggplot2", "ggthemes", "jtools", "tidyverse", "scales")

```

#### 20 year group

```{r}

Panel_long <- read_xlsx("//tsclient/C/Users/tanay/Downloads/Sims/Panel_28Jan2024.xlsx")
Panel_long <- read_xlsx("C:/Users/tanay/Downloads/Sims/Panel_28Jan2024.xlsx")

Avg_1_20yrs <- as.numeric(Panel_long %>% filter(ADS_ever == 0) %>% summarise(Avg1 = mean(EDI)))   
Avg_2_20yrs <- as.numeric(Panel_long %>% filter(ADS_ever == 1, ADS2 == 0, ADS_IN_0_1_2 == 0) %>% summarise(Avg2 = mean(EDI)))    
Avg_3_20yrs <- as.numeric(Panel_long %>% filter(ADS_ever == 1, ADS2 == 0, ADS_IN_0_1_2 == 1) %>% summarise(Avg3 = mean(EDI)))

```


```{r}
library(knitr)

# Render each table as HTML using kable()
tab1_html <- kable(Avg_1_20yrs, caption = "ADS_ever = 0")
tab2_html <- kable(Avg_2_20yrs, caption = "ADS_ever = 1, ADS_IN_0_1_2 = 0, ADS2 = 0")
tab3_html <- kable(Avg_3_20yrs, caption = "ADS_ever = 1, ADS_IN_0_1_2 = 1, ADS2 = 0")

tab1_html
tab2_html
tab3_html
```


```{r}
df <- select(Panel_long, code, year, ADS2, ADS_ever,  ADS_IN_0_1_2, EDI)

df <- mutate(df, EDI_group_avg = case_when(
  ADS_ever == 0 ~ Avg_1_20yrs,
  ADS_ever == 1 & ADS2 == 0 & ADS_IN_0_1_2 == 0 ~ Avg_2_20yrs,
  ADS_ever == 1 & ADS2 == 0 & ADS_IN_0_1_2 == 1 ~ Avg_3_20yrs,
  TRUE ~ NA_real_
))  %>%
  mutate(EDI_group_avg = replace_na(EDI_group_avg, 0)) %>% select(code, year, EDI_group_avg, EDI)

```


```{r}


simulate_poisson_count <- function(lambda, window_end = 1) {
  if (is.na(lambda) || lambda <= 0) {
    return(0)
  }

  time <- 0
  count <- 0

  while (time < window_end) {
    time <- time + rexp(1, rate = lambda)
    if (time < window_end) {
      count <- count + 1
    }
  }

  return(count)
}


Simulated_regs <- list()
U <- c(0.02, 0.05, 0.1, 0.25, 0.5, 0.75)


for(u in U){
  ##    df <- mutate(results_df, lambda = Pred/(1-u))
  df <- mutate(df, lambda = EDI_group_avg*u)
  
  # preallocate a list of length 1000 for this u
  Simulated_regs[[as.character(u)]] <- vector("list", 1000)
  
  for(i in 1:1000){
    df_sim <- df %>%  rowwise() %>%
      mutate(
        EDI_new_scaled = 100 * (EDI + simulate_poisson_count(lambda, 1))
      ) %>% ungroup()
    
    df_sim <- left_join(Panel_long, df_sim, by = c("code", "year")) 
    
    reg <- feols(
      EDI_new_scaled ~ ADS2 + ADS_IN_0_1_2 + neigh + HD_dummy_0_1 + dist | code + year,
      data = df_sim,
      conley(25)
    )
    
    # save coefficient + se as a pair
    Simulated_regs[[as.character(u)]][[i]] <- c(
      coef = reg$coefficients[1],
      se   = reg$se[1]
    )
  }
}

```

```{r}
dat <- do.call(rbind, lapply(names(Simulated_regs), function(cat) {
  tmp <- do.call(rbind, Simulated_regs[[cat]])  # bind 1000 x 3 matrix
  data.frame(
    category = cat,
    coefficient = tmp[, 1],
    se = tmp[, 2]
  )
}))

dat <- mutate(dat, UR_perc = as.factor(round(as.numeric(category)/(1+ as.numeric(category)), 2)))

# Compute mean coefficient and empirical SE
dat_summary <- dat %>%
  group_by(UR_perc) %>%
  summarise(
    mean_coef = mean(coefficient),
    emp_se = sd(coefficient) / sqrt(n()),
    .groups = "drop"
  ) %>%
  mutate(
    lower = mean_coef - 1.96 * emp_se,
    upper = mean_coef + 1.96 * emp_se
  )

```

```{r}
dat_summary

# Plot
p <- ggplot(dat, aes(x = UR_perc, y = coefficient, fill = UR_perc)) +
  geom_boxplot(width = 0.1, outlier.size = 0.5, alpha = 0.7) +
  theme_clean(base_size = 18) +
  theme(legend.position = "none") +
  labs(
    title = "Simulated effects under undderreported\n elephant deaths",
    x = "% deaths underreported",
    y = "Coefficient"
  )

p
```



```{r}
Panel_long <- Panel_long %>%
  group_by(code) %>%  mutate(ADS_stag = cummax(ADS2)) %>%
  ungroup()



Reg_stag <- list(

feols(EDI_scaled ~ ADS_stag + ADS_IN_0_1_2  + neigh +  HD_dummy_0_1 + dist | code + year, data = Panel_long, conley(25)),

feols(EDI_accident_scaled ~ ADS_stag + ADS_IN_0_1_2  + neigh +  HD_dummy_0_1 + dist | code + year, data = Panel_long, conley(25)),

feols(EDI_other_scaled ~ ADS_stag + ADS_IN_0_1_2  + neigh +  HD_dummy_0_1 + dist | code + year, data = Panel_long, conley(25)))

```

```{r}

##Table: EDI
etable(Reg_stag, title = "Effects of Anti-Depredation Squads on elephant deaths (Staggered treatments)", fitstat = c("n", "aic", "bic"), style.tex = style.tex(main = "base"), replace = T, signif.code = c("***"=0.01, "**"=0.05, "*"=0.10), drop.section = c("fixef"), extralines = list("Fixed effects" = c("Two-way", "Two-way", "Two-way"), "Standard Errors" = c("Conley", "Conley", "Conley")), dict = c("Pop" = "Population", "ADS_stag" = "ADS_stag", "dist" = "dist_eleph_habitat", "neigh" = "ADS_neighbor", "HDI_scaled" = "HD", "EDI.wo.Hits_scaled" = "ED w/o hits", "EDI.wo.Hits" = "ED w/o hits", "HD_dummy_0_1" = "HD_0_1", "EDI_accident_scaled" = "EDI (accident)", "EDI_accident" = "EDI (accident)", "EDI_other_scaled" = "EDI (other)", "EDI_other" = "EDI (other)"))

```
